# ------------------------------------------------------------------------------#
# Description: Create Table A.7: Robustness of peer effects to spatial fixed effects
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(data.table)
library(fixest)
library(modelsummary)

options(scipen = 999)

# Load adoption data -----------------------------------------------------------
load("data/adoptions/adoptions.RData")

# Merge peer adoption buckets --------------------------------------------------
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

# Merge peer eligibility buckets -----------------------------------------------
load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

# Merge total peers ------------------------------------------------------------
load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

# Merge parcel data ------------------------------------------------------------
load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

# Merge census data ------------------------------------------------------------
load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale     := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale  := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale      := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale      := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale    := peer.adopt.both.100 / 100]

# Keep post-2010 only ----------------------------------------------------------
dat <- dat[year > 2010]

# Create average eligibility variable ------------------------------------------
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Merge shoreline data ---------------------------------------------------------
load("data/shoreline/shoreline.RData")
setnames(near.shore, "near_dist", "near_shore")
dat <- merge(dat, near.shore[, .(pin, near_shore)], by = "pin", all.x = TRUE)
rm(near.shore)

# Merge arterial road data -----------------------------------------------------
load("data/arterial/arterial.RData")
setnames(near.arterial, "near_dist", "near_road")
dat <- merge(dat, near.arterial[, .(pin, near_road)], by = "pin", all.x = TRUE)
rm(near.arterial)

# Merge mandatory GSI data -----------------------------------------------------
load("data/mandatory_gsi/peer_adopt_buckets_any_mandgsi.RData")
dat <- merge(dat, peer.adopt.bucket.any.mandgsi, by = c("pin", "year"))
rm(peer.adopt.bucket.any.mandgsi)

# Merge DEM data ---------------------------------------------------------------
load("data/dem/all_parcels_dem.RData")
dat <- merge(dat, parcels.dem, by = "pin")
rm(parcels.dem)

# Estimate models --------------------------------------------------------------
m1 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              blkgrp + year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m2 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              basin + year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "basin", data = dat[elig.cs == 1 & post.rainwise == 0])

m3 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              zip5 + year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "zip5", data = dat[elig.cs == 1 & post.rainwise == 0])

m4 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              sub.area + year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "sub.area", data = dat[elig.cs == 1 & post.rainwise == 0])

m5 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              tract + year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "tract", data = dat[elig.cs == 1 & post.rainwise == 0])

m6 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              basin^year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "basin", data = dat[elig.cs == 1 & post.rainwise == 0])

m7 <- feols(adopt.any * 100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob | 
              blkgrp^year | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

# Load IV diagnostics ----------------------------------------------------------
load("data/iv/iv_diag_spatial.RData")
for (i in 1:7) {
  assign(paste0("eff.f.", i), round(iv_diag_spatial[[i]]$F_stat[4], 2))
  assign(paste0("ar.ci.", i), paste0("[", round(iv_diag_spatial[[i]]$AR$ci[1], 2), ", ",
                                     round(iv_diag_spatial[[i]]$AR$ci[2], 2), "]"))
}

# Define Variable Names ----------------------------------------------------------
varnames <- c('year' ='Year', 'area' = 'Area','zip5' = 'Zipcode','basin'='Basin',
              'sub.area' = 'Sub-Area', 'tract' = 'Tract','blkgrp' = 'Block Group',
              'blk' = 'Block', 'peer.adopt.any.100'= '$\\widehat{\\text{Peer Adoptions}}$')


# Preview Table A.7 ------------------------------------------------------------
style_tex <- style.tex(tablefoot.value = '',
                       model.title = '',
                       depvar.title = '',
                       fixef.title = '\\midrule \\emph{Fixed effects, controls, and sample}',
                       var.title = '\\midrule',
                       stats.title = "\\midrule")

etable(m1, m2, m3, m4, m5, m6, m7,
       keep = c("%peer.adopt.any.100"),
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       fixef_sizes.simplify = FALSE, fitstat = c("n"),
       extralines = list(
         "_F" = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5, eff.f.6, eff.f.7),
         "_AR CI" = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5, ar.ci.6, ar.ci.7)
       ))

# Save Table A.7 ---------------------------------------------------------------
etable(m1, m2, m3, m4, m5, m6, m7,
       keep = c("%peer.adopt.any.100"),
       tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       fixef_sizes.simplify = TRUE,
       fitstat = c("n", "ivfall"), digits = 2,
       style.tex = style_tex,
       extralines = list(
         "_F" = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5, eff.f.6, eff.f.7),
         "_AR CI" = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5, ar.ci.6, ar.ci.7)
       ),
       file = "output/tables/table_a7.tex", replace = TRUE)
